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[i]  Estimation  of  the  three-dimensional  (3-D)  distribution  of  hydrologic  properties 
and  related  uncertainty  is  a  key  for  improved  predictions  of  hydrologic  processes  in  the 
subsurface.  However  it  is  difficult  to  gain  high-quality  and  high-density  hydrologic 
information  from  the  subsurface.  In  this  regard  a  promising  strategy  is  to  use  high- 
resolution  geophysical  data  (that  are  relatively  sensitive  to  variations  of  a  hydrologic 
parameter  of  interest)  to  supplement  direct  hydrologic  information  from  measurements  in 
wells  (e.g.,  logs,  vertical  profiles)  and  then  generate  stochastic  simulations  of  the 
distribution  of  the  hydrologic  property  conditioned  on  the  hydrologic  and  geophysical  data. 

In  this  study  we  develop  and  apply  this  strategy  for  a  3-D  field  experiment  in  the 
heterogeneous  aquifer  at  the  Boise  Hydrogeophysical  Research  Site  and  we  evaluate  how 
much  benefit  the  geophysical  data  provide.  We  run  high-resolution  3-D  conditional 
simulations  of  porosity  with  both  simulated-annealing-based  and  Bayesian  sequential 
approaches  using  information  from  multiple  intersecting  crosshole  gound-penetrating  radar 
(GPR)  velocity  tomograms  and  neutron  porosity  logs.  The  benefit  of  using  GPR  data  is 
assessed  by  investigating  their  ability,  when  included  in  conditional  simulation,  to  predict 
porosity  log  data  withheld  from  the  simulation.  Results  show  that  the  use  of  crosshole  GPR 
data  can  significantly  improve  the  estimation  of  porosity  spatial  distribution  and  reduce 
associated  uncertainty  compared  to  using  only  well  log  measurements  for  the  estimation. 

The  amount  of  benefit  depends  primarily  on  the  strength  of  the  petrophysical  relation 
between  the  GPR  and  porosity  data,  the  variability  of  this  relation  throughout  the 
investigated  site,  and  lateral  structural  continuity  at  the  site. 

Citation:  Dafflon,  B.,  and  W.  Barrash  (2012),  Three-dimensional  stochastic  estimation  of  porosity  distribution:  Benefits  of  using 
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1.  Introduction 

[2]  Detailed  knowledge  of  the  distribution  of  parameters 
controlling  flow  and  transport  processes  in  aquifers  is  a  key 
prerequisite  for  realistic  simulation  of  groundwater  flow 
and  contaminant  transport,  and  thus  for  sustainable  man¬ 
agement  and  effective  remediation  of  groundwater  resour¬ 
ces.  Estimation  of  aquifer  parameters  is  difficult,  however, 
given  the  challenges  with  direct  sampling  of  the  subsurface, 
and  the  practical  limitations  on  capturing  three-dimensional 
(3-D)  spatial  structure.  In  this  regard  one  promising 
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strategy  is  to  use  high-resolution  geophysical  data  (that  are 
relatively  sensitive  to  variations  of  a  hydrologic  parameter 
of  interest)  to  supplement  direct  hydrologic  information 
from  measurements  in  wells  (e.g.,  logs,  vertical  profiles) 
and  then  generate  stochastic  simulations  of  the  distribution 
of  the  hydrologic  property  conditioned  on  the  hydrologic 
and  geophysical  data.  Geophysical  data  can  provide  infor¬ 
mation  at  spatial  scales  and  locations  that  are  unattainable 
with  conventional  hydrologic  measurement  techniques, 
given  a  functional  or  statistical  relation  between  geophysi¬ 
cal  and  hydrologic  parameters  [e.g.,  Hubbard  and  Rubin, 
2005].  In  this  way  geophysical  data  have  successfully 
assisted  hydrologic  investigations  with  (1)  mapping  of 
hydrogeologically  relevant  structures  [e.g.,  Beres  et  al., 
1995;  Asprion  and  Aigner,  1999;  Barrash  and  Clemo, 
2002],  (2)  hydrologic  parameter  estimation  [e.g.,  Hyndman 
et  al.,  2000;  Chen  et  al.,  2001 ;  Kowalskv  et  al,  2005 ;  Linde 
et  al.,  2006;  Paasche  et  al.,  2006;  Harp  et  al.,  2008; 
Dafflon  et  al.,  2010;  Straface  et  al.,  2011],  and  (3)  visualiza¬ 
tion  of  subsurface  processes  by  monitoring  temporal  hydro- 
logic  changes  in  the  subsurface  [e.g.,  Kemna  et  al,  2002; 
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Day-Lewis  et  al.,  2003 ;  Singha  and  Gorelick,  2005 ;  Johnson 
el  al.,  2007], 

[3]  Integration  of  geophysical  and  hydrologic  data  with 
stochastic  simulation  approaches,  such  as  sequential  Gaus¬ 
sian  and  Bayesian  approaches  and  Monte-Carlo  (MC)-type 
methods  [e.g.,  Deutsch,  2002;  Kelkar  and  Perez,  2002], 
produce  multiple  realizations  of  a  target  parameter  that  are 
consistent  with  available  measurements  and  spatial  structure 
information.  An  advantage  of  such  approaches  is  that  the  set 
of  simulated  realizations  represents  the  joint  uncertainty  in 
our  knowledge  of  subsurface  properties.  For  example  se¬ 
quential  Gaussian  simulation  [e.g.,  Deutsch  and  Journel, 
1998]  has  been  used  successfully  to  integrate  crosshole  seis¬ 
mic  and  geostatistical  information  by  Hyndman  et  al.  [2000], 
and  also  by  McKenna  and  Poeter  [1995]  for  generating  indi¬ 
cator  simulations  based  on  crosshole  seismic  tomography 
and  borehole  logging  information.  Bayesian  approaches 
[e.g.,  Gelman  et  al.,  2003]  have  been  applied  to  integrate 
diverse  geophysical  and  hydrologic  data  [e.g.,  Ezzedine 
et  al.,  1999;  Chen  et  al.,  2001],  and  can  be  included  within  a 
sequential  simulation  approach  [e.g.,  Scheibe  and  Chien, 
2003].  Recently  Dubreuil-Boisclair  et  al.  [2011]  successfully 
used  Bayesian  sequential  simulation  to  integrate  crosshole 
GPR  velocity  and  amplitude  tomograms  with  well  measure¬ 
ments  to  image  the  hydraulic  conductivity  distribution. 
Also,  another  simulation  strategy  that  is  particularly  flexible 
for  integrating  various  data  sets,  constraints,  and  petrophysi¬ 
cal  relations  is  the  optimization  of  a  parameter  realization 
through  a  Monte-Carlo  (MC)  perturbation  process  for  a  given 
objective  function  [e.g.,  Deutsch  and  Wen,  2000;  Kelkar  and 
Perez,  2002;  Tronicke  and  Holliger,  2005;  Dafflon  et  al., 
2009a],  Among  the  MC  methods,  a  simulated-annealing 
(SA)-based  approach  developed  by  Dafflon  et  al.  [2009a] 
generates  realistic  spatial  distributions  of  porosity  using 
crosshole  GPR  tomograms  and  porosity  log  data. 

[4]  The  objective  of  this  study  is  to  develop  improved 
stochastic  integration  techniques  for  3-D  estimation  of  the 
porosity  distribution,  and  to  demonstrate  the  approaches 
with  data  from  the  Boise  Hydrogeophysical  Research  Site 
(BHRS).  We  evaluate  added  benefit  from  using  multidirec¬ 
tional  crosshole  GPR  velocity  tomograms  [ Dafflon  et  al., 
2011a]  with  porosity  log  data  [Barrash  and  Clemo,  2002] 
for  the  estimation  of  the  3-D  distribution  of  porosity  and 
associated  uncertainty.  The  basis  to  relate  GPR  velocity  to 
porosity  is  the  strong  sensitivity  of  GPR  velocity  to  soil 
water  content,  which  is  equivalent  to  porosity  in  the  satu¬ 
rated  zone  [e.g.,  Schon,  2004],  High-resolution  conditional 
simulations  of  the  porosity  spatial  distribution  using  these 
data  are  generated  using  both  SA-based  [e.g.,  Dafflon  et  al., 
2009a]  and  Bayesian  simulation  approaches  [e.g.,  Dubreuil- 
Boisclair  et  al.,  2011].  We  selected  these  methods  because 
of  their  demonstrated  success  compared  to  other  approaches 
in  2-D  conditional  simulations  of  hydrologic  parameter 
distributions  using  data  similar  to  those  used  in  this  study 
[ Dafflon  et  al.,  2009b;  Dubreuil-Boisclair  et  al.,  2011].  Here 
these  approaches  are  adapted  to  3-D  and  to  the  general  case 
of  variably  spaced  data  (as  are  available  at  the  BHRS  and 
many  sites).  The  benefit  added  by  including  crosshole  GPR 
tomograms  in  the  stochastic  estimation  of  the  porosity  spa¬ 
tial  distribution  is  assessed  by  evaluating  how  well  porosity 
values  are  predicted  at  locations  where  porosity  log  data  are 
withheld.  This  is  done  for  two  cases  where  the  withheld 


porosity  log  is  located:  (1)  at  a  crosshole  GPR  profile,  so  the 
full  available  data  set,  with  the  exception  only  of  the  with¬ 
held  porosity  log,  is  used  to  condition  the  3-D  simulation; 
and  (2)  where  no  data  are  available  so  all  the  existing  cross¬ 
hole  profiles  intersecting  the  withheld  porosity  log  are  also 
withheld  for  this  estimation  process.  For  both  cases  we  do 
this  for  two  wells  having  locally  variable  stratigraphy. 

[5]  Developing  an  effective  integration  approach  to  esti¬ 
mate  the  porosity  spatial  distribution  in  unconsolidated 
aquifers  is  important  because  porosity  heterogeneity  can 
play  a  significant  role  in  transport  behavior  [Hassan,  2001 ; 
Hu  et  al.,  2009;  Dafflon  et  al,  2011b;  Li  et  al.,  2012]. 
Also,  although  the  spatial  distribution  of  hydraulic  conduc¬ 
tivity  generally  is  much  more  difficult  to  estimate  than  po¬ 
rosity,  these  parameters  commonly  exhibit  some  degree  of 
similarity  with  regard  to  their  spatial  variability  and/or  spa¬ 
tial  correlation  [e.g.,  Hubbard  et  al.,  2001;  Scheibe  and 
Chien,  2003]. 

[6]  In  this  paper  we  first  describe  the  simulation  approaches 
used  to  perform  3-D  conditional  simulations.  Next,  we  briefly 
describe  the  porosity  log  data  and  GPR  velocity  tomograms 
from  the  BHRS,  and  the  conditional  petrophysical  relations 
and  geostatistical  models  used  in  the  simulation  process.  Then 
we  evaluate  the  benefit  gained  from  including  crosshole  GPR 
data  in  the  simulation  process  by  predicting  porosity  with  sev¬ 
eral  approaches.  Finally,  we  discuss  improvements  achieved 
in  imaging  hydrostratigraphic  units  at  the  BHRS. 

2.  Stochastic  Approaches 

2.1.  SA-based  Optimization  Approach 

[7]  SA  is  a  directional  MC-type  optimization  procedure 
[Deutsch  and  Wen,  1998;  Day-Lewis  et  al.,  2000;  Kelkar 
and  Perez,  2002;  Tronicke  and  Holliger,  2005;  Dafflon 
et  al.,  2009a]  that  involves  the  repeated  perturbation  of  a 
random  field  in  order  to  satisfy  a  multicomponent  global 
objective  function.  In  this  study  we  use  the  SA-based  condi¬ 
tional  simulation  approach  of  Dafflon  et  al.  [2009a]  which 
was  developed  to  simulate  hydrologic  parameter  fields  by 
assimilating  both  larger-scale  subsurface  structures  (pro¬ 
vided  by  geophysical  data)  and  smaller-scale  fluctuations 
(provided  by  well  log  data  and  spatial  correlation  functions). 
This  approach  consists  of:  (1)  generation  of  an  initial  real¬ 
ization  of  the  target  parameter  (e.g.,  porosity),  (2)  random 
selection  and  perturbation  of  a  cell  in  the  model  by  drawing 
from  a  probability  distribution  of  the  target  parameter  given 
the  available  data  at  that  location  (e.g.,  using  conditional  pet¬ 
rophysical  distributions),  (3)  updating  the  objective  function 
(e.g.,  spatial  correlation  functions),  (4)  accepting  the  pertur¬ 
bation  if  the  objective  function  is  lowered  according  to  a 
Boltzmann-type  exponential  probability  function,  (5)  repeti¬ 
tion  of  steps  (2)  to  (4)  until  the  objective  function  is  satisfied. 
The  initial  realization  is  generated  by  performing  the  second 
step  for  each  cell  in  the  model  to  ensure  that  each  initial  real¬ 
ization  is  different. 

[8]  Here  we  adapted  the  approach  of  Dafflon  et  al. 
[2009a]  to  3-D  and  to  variably  spaced  data.  To  this  end,  the 
generation  of  the  initial  realization  (1)  and  each  perturba¬ 
tion  (2)  depends  on  the  available  data  at  the  considered 
cell.  That  is,  if  crosshole  GPR  data  are  available  at  a  given 
cell,  the  new  value  is  drawn  randomly  from  a  probability 
distribution  of  porosity  given  the  available  GPR  data  at  the 
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cell  location  using  the  conditional  petrophysical  distribu¬ 
tion  developed  from  collocated  GPR  and  porosity  data  at 
the  wells  (also  see  discussion  in  section  4.1).  Next,  at  and 
near  the  boreholes  the  probability  distribution  of  the  poros¬ 
ity  is  more  tightly  constrained  with  a  conditional  expecta¬ 
tion  equal  to  the  porosity  log  data  and  a  relatively  small 
standard  deviation  defining  its  uncertainty.  Finally,  at  cells 
where  no  data  are  available  the  conditional  expectation  is 
obtained  through  kriging  the  log  data  and  the  conditional 
expected  porosity  values  inferred  from  the  crosshole  GPR 
data  through  the  conditional  petrophysical  distribution.  The 
variance  of  the  porosity  probability  distribution  at  such 
cells  is  set  equal  to  the  variance  of  all  the  log  data  present 
in  a  depth  interval  centered  on  the  simulated  location  and 
equal  to  the  range  of  the  vertical  spatial  correlation  func¬ 
tion  of  the  porosity.  We  do  this  instead  of  using  the  kriging 
variance  because  it  is  difficult  to  assess  the  true  variance  at 
these  locations,  and  so  we  prefer  to  allow  relatively  large 
variability  where  no  data  are  available  (i.e.,  to  overestimate 
versus  underestimate  the  true  uncertainty). 

[9]  The  objective  function  (3)  is  defined  by  matching  a 
parametric  spatial  correlation  function  in  x,  y  and  z  direc¬ 
tions.  Dafflon  et  al.  [2009a]  have  shown  that  these  correla¬ 
tion  models  need  to  be  considered  only  at  “relatively 
small”  lags  where  the  GPR  data  fail  to  capture  the  hetero¬ 
geneity  and  where  variability  can  be  estimated  relatively 
well  from  available  porosity  log  data,  while  at  large  scales 
the  variability  is  constrained  only  by  GPR  data.  This  is  sim¬ 
ilar  to  assigning  a  very  high  uncertainty  to  the  porosity  spatial 
correlation  function  at  large  lags.  Finally,  a  last  parameter 
that  can  be  important  in  the  optimization  process  is  the  “tem¬ 
perature”  parameter  in  the  Boltzmann-type  exponential  prob¬ 
ability  function  [e.g.,  Dafflon  et  al.,  2009a],  which  controls 
the  extent  to  which  perturbations  that  do  not  lower  the  objec¬ 
tive  function  are  accepted  (4).  In  this  study,  due  to  the  rela¬ 
tive  simplicity  of  the  objective  function  and  based  on  some 
preliminary  tests,  the  temperature  parameter  does  not  signifi¬ 
cantly  influence  the  process  and  therefore  has  been  set  rela¬ 
tively  low. 

2.2.  Bayesian  Sequential  Simulation 

[10]  The  Bayesian  sequential  simulation  approach  is 
mainly  based  on  traditional  sequential  Gaussian  simulation 
[e.g.,  Deutsch  and  Journel,  1998]  with  the  added  use  of 
Bayesian  formula  [e.g.,  Chen  et  al.,  2001;  Gelman  et  al., 
2003].  The  Bayesian  methodology  has  been  used  widely 
with  success  for  updating  prior  knowledge  with  additional 
conditional  sources  of  information  [e.g.,  Ezzedine  et  al., 
1999;  Chen  et  al.,  2001 ;  Linde  et  al.,  2007],  and  especially 
in  the  context  of  Bayesian  sequential  simulations  [Gastaldi 
et  al.,  1998;  Scheibe  and  Chien,  2003;  Dubreuil-Boisclair 
et  al.,  2011].  In  this  approach  each  new  simulated  value  is 
drawn  from  a  probability  distribution  obtained  by  updating 
a  prior  distribution  of  a  target  parameter  into  a  posterior 
distribution  given  additional  data  at  a  given  location  and 
using  estimates  of  the  conditional  petrophysical  relation 
between  collocated  geophysical  and  hydrologic  data.  In  a 
recent  study  by  Dubreuil-Boisclair  et  al.  [2011]  the  2-D  hy¬ 
draulic  conductivity  distribution  between  two  boreholes 
was  estimated  with  Bayesian  sequential  simulation  using 
collocated  results  from  a  GPR  velocity  and  attenuation 
tomogram  and  hydraulic  conductivity  measurements  at  the 


wells  to  estimate  multivariable  relations  between  hydraulic 
conductivity,  electrical  conductivity,  and  permittivity  with 
a  kernel  probability  function  [e.g.,  Wand  and  Jones,  1995]. 
Their  Bayesian  simulation  approach  showed  clear  improve¬ 
ment  compared  to  other  approaches  such  as  a  sequential 
co-simulation  approach. 

[11]  The  3-D  Bayesian  sequential  simulation  approach 
used  here  consists  of  two  stages.  The  first  stage  is  per¬ 
formed  only  at  cells  located  along  crosshole  GPR  profiles. 
The  process  consists  of  (1)  selecting  a  cell  that  has  not  been 
visited  yet,  (2)  applying  simple  kriging  under  a  Gaussian 
assumption  using  all  the  available  values  of  porosity  to  esti¬ 
mate  the  expected  mean  and  variance  that  define  the  prior 
distribution  at  the  visited  cell,  (3)  defining  the  likelihood 
distribution  obtained  from  both  the  available  crosshole 
GPR  data  at  the  selected  cell  and  the  conditional  petrophys¬ 
ical  distribution  inferred  from  the  collocated  crosshole 
GPR  and  porosity  log  data  at  the  wells,  (4)  inferring  the 
posterior  probability  distribution  by  applying  Bayes’  rule 
to  update  the  prior  distribution  with  the  likelihood  distribu¬ 
tion  (by  calculating  their  product),  (5)  drawing  a  porosity 
value  from  the  posterior  probability  distribution  for  the 
selected  cell  and  adding  it  to  the  available  porosity  data 
before  visiting  the  next  cell,  and  (6)  repeating  the  above 
steps  until  all  the  selected  cells  have  been  visited.  It  is 
worth  noting  that  the  conditional  petrophysical  distribution 
used  to  estimate  the  likelihood  distribution  (3)  is  identical 
to  the  one  used  in  the  SA  approach  (see  also  section  4.1). 

[12]  The  second  step  of  our  approach  is  to  perform  stand¬ 
ard  sequential  Gaussian  simulation  in  all  the  remaining 
cells  that  do  not  contain  any  data.  In  fact  this  stage  is  iden¬ 
tical  to  the  previous  one  except  that  no  likelihood  function 
is  available  and  thus  the  prior  distribution  alone  defines  the 
posterior  probability  distribution  from  which  a  value  of  po¬ 
rosity  is  drawn.  Finally  we  note  that,  although  we  use  sim¬ 
ple  kriging  in  this  paper,  we  have  conducted  additional 
investigations  that  show  the  use  of  ordinary  kriging  does 
not  change  the  results  of  this  study  significantly. 

3.  Hydrogeologic  Field  Site 

3.1.  Well  Measurements 

[13]  The  BHRS  is  a  hydrologic  and  geophysical  field 
research  site  located  near  Boise,  Idaho,  USA.  The  shallow 
unconfined  aquifer  at  the  site  resides  within  an  approxi¬ 
mately  20-m-thick  layer  of  sediments  consisting  of  coarse, 
unconsolidated,  fluvial  deposits  [Barrash  and  Clemo, 
2002 ;  Barrash  and  Reboulet,  2004]  with  minimal  fractions 
of  silt  and  clay,  and  is  underlain  by  a  layer  of  red  clay.  At 
this  site  1 8  wells  have  been  emplaced  in  a  configuration  to 
facilitate  a  wide-variety  of  hydrologic  and  geophysical  test¬ 
ing  [. Barrash  et  al.,  1999].  All  of  the  wells  were  cased  with 
4-in  PVC  well  screen  and  were  carefully  completed  to  mini¬ 
mize  the  disturbance  of  the  surrounding  formation.  The  well 
field  consists  of  13  wells  in  the  central  area  (~20  m  in  diam¬ 
eter)  and  five  boundary  wells  about  10  to  35  m  from  the  cen¬ 
tral  area  (Figure  1).  The  general  design  of  the  13  inner  wells 
is  a  central  well  (Al)  surrounded  by  two  concentric  rings  of 
six  wells  each  (B1-B6  and  C1-C6,  respectively).  The  distan¬ 
ces  between  these  wells  vary  between  2.6  and  8.6  m. 

[14]  Porosity  values  in  each  borehole  were  estimated  at 
0.06  m  intervals  from  the  measured  neutron-log  count  rate 
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Figure  1.  Aerial  view  of  the  BHRS  and  detailed  map  of 
the  central  area  wells  with  lines  locating  crosshole  GPR 
tomograms.  Joint  inversion  of  multiple  GPR  profiles  to 
reconstruct  a  multidirectional  model  of  GPR  velocity  was 
performed  for  both  31  (continuous  gray  and  black  lines) 
and  seven  data  sets  together  (black  lines). 

[Barrash  and  Clemo,  2002]  through  a  petrophysical  trans¬ 
form  [Hearst  and  Nelson,  1985]  calibrated  with  measure¬ 
ments  of  porosity  in  similar  environments  [ Barrash  and 
Clemo,  2002].  Based  on  these  porosity  logs  Barrash  and 
Clemo  [2002]  identified  five  units  with  distinct  spatial  dis¬ 
tribution,  mean,  and  variance.  The  fitted  Gaussian  functions 
to  the  porosity  distributions  of  Units  1  to  5  have  means 
equal  to  0.18,  0.24,  0.172,  0.224  and  0.425,  respectively, 
and  standard  deviations  equal  to  0.022,  0.038,  0.024,  0.05 
and  0.055,  respectively.  Unit  5  is  a  high  porosity  sand  chan¬ 
nel  that  thickens  toward  the  Boise  River  and  pinches  out  in 
the  center  of  the  well  field.  Units  1-4  are  a  sequence  of  con¬ 
glomerates  with  gravel  and  cobble  framework  and  sand  to 
pebble  matrix  in  the  interstices  of  the  framework.  More 
recently  capacitive  electrical  conductivity  measurements 
have  identified  a  subunit  (Unit  2b)  that  is  present  in  all  the 
wells  shown  in  the  inset  of  Figure  1  except  wells  Bl,  B3, 
Cl  and  C2  [Mwenifunibo  et  al.,  2009],  Four  of  the  five 
main  units  (i.e.,  Units  1-4)  occur  in  the  depth  interval  dis¬ 
cussed  in  this  paper  (Figure  2). 

3.2.  Crosshole  GPR  Tomograms 

[is]  A  total  of  38  crosshole  GPR  data  sets  were  acquired 
from  1998  to  2000  between  the  13  wells  in  the  central  area 
of  the  BFIRS  (Figure  1)  using  a  Mala  Ramac  GPR  system 
with  antennas  having  a  nominal  center  frequency  in  air  of 
250  MFlz.  For  this  study  we  used  31  crosshole  data  sets 


Figure  2.  Velocity  tomogram  for  crosshole  profile 
between  A1-B2  from  inversion  of  seven  crosshole  GPR 
data  sets  together  (modified  from  Dafflon  et  al.  [2011a]). 
Intersections  of  other  profiles  through  the  A1-B2  profile  are 
delineated  (gray).  The  correlation  between  velocity  and 
porosity-log  data  at  the  borehole  location  is  given  at  the 
bottom  of  each  log  comparison  panel.  For  comparison  with 
stratigraphy  at  the  BFIRS,  contacts  between  Units  1  to  4  are 
identified  in  the  porosity  logs  [. Barrash  and  Clemo,  2002 ; 
Mwenifunibo  et  al.,  2009]. 

between  1 1  wells.  All  the  data  have  a  common  depth¬ 
sampling  interval  of  approximately  0.2  m  in  the  transmitter 
and  receiver  boreholes  and  cover  the  saturated  zone 
between  846.5  and  832  m  elevation.  Antenna  positions 
were  corrected  to  account  for  borehole  deviations  using 
measurements  from  a  magnetic  deviation  logging  survey  in 
2010.  Travel  times  of  the  direct  transmitted  wavefield  for 
each  crosshole  GPR  data  set  were  determined  using  a  semi- 
automated  picking  procedure. 

[16]  GPR  velocity  tomograms  used  in  this  study  were 
obtained  by  jointly  inverting  intersecting  crosshole  GPR 
travel  time  data  sets  [. Dafflon  et  al.,  2011a].  In  this  method, 
first  all  transmitter  and  receiver  coordinates,  and  coordi¬ 
nates  of  the  related  profile  intersections,  are  projected  in 
2-D.  Then  the  tomographic  kernel  matrix,  regularization 
matrix,  and  model  vector  are  constructed  such  that  ( 1 )  com¬ 
mon  cells  are  considered  where  profiles  intersect  and  (2)  all 
profiles  are  considered  together  in  a  single  expanded  tomo¬ 
graphic  system.  The  velocity  is  inverted  for  the  set  of  all 
2-D  cells  through  ray-based  inversion  using  a  conjugate- 
gradient-based  least-squares  approach  [e.g.,  Scales,  1987] 
and  accounting  for  possible  errors  in  the  crosshole  GPR 
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survey  geometry,  and  for  antenna  transmission  behavior 
and  locations  [Irving  et  al,  2007 ;  Dcifjlon  et  al.,  201  la]. 

[17]  The  above  inversion  strategy  has  been  shown  to  be 
a  robust  procedure  for  two  different  cases  at  the  BHRS: 
(1)  31  profiles,  or  all  the  available  data  with  the  exception  of 
the  profiles  obtained  using  wells  C4  and  C5,  and  (2)  a  subset 
of  seven  profiles  of  higher-than-average  quality  in  a  corre¬ 
spondingly  limited  area  of  the  whole  data  set  (Figure  1). 
The  GPR  velocity  tomograms  inverting  seven  profiles 
together  (from  Dafflon  et  al.  [2011a])  and  shown  here  for 
A1-B2  (Figure  2):  (1)  successfully  image  the  major  strati¬ 
graphic  units  and  their  contacts  at  the  BHRS,  and  (2)  con¬ 
firm  the  variability  differences  within  and  between  units 
evident  in  porosity  data  [Barrash  and  Clemo,  2002], 
including  recognition  of  smaller-scale  subfacies  structures. 
The  correlation  coefficients  between  the  collocated  GPR 
velocities  and  porosity  data  at  the  11  wells  used  in  this 
study  (Figure  1)  range  between  —0.32  and  —0.79  with  a 
mean  of  —0.57.  We  note  that  there  is  no  link  between  the 
values  of  the  correlation  coefficient  and  the  number  of 
GPR  crosshole  profiles  intersecting  at  each  well  location 
[ Dafflon  et  al.,  2011a].  With  few  exceptions,  these  correla¬ 
tions  can  be  considered  as  significant  given  that  GPR  data 
have  a  more  limited  spatial  resolution  than,  and  are  not  able 
to  define  small-scale  heterogeneity  near  boreholes  as  well 
as,  the  porosity  measurements  (Figure  2).  Also  these  corre¬ 
lation  coefficients  are  known  to  be  relatively  decreased  by 
the  presence  of  subunit  2b  which  is  less  well  imaged  by  the 
GPR  data  due  to  anomalous  conductivity/dielectric  petro¬ 
physical  behavior  not  dominated  by  water-saturated  poros¬ 
ity  [Mwenifumbo  et  ah,  2009 ;  Dafflon  et  al.,  2011a]. 

4.  Conditional  Petrophysical  Distribution  and 
Spatial  Correlation  Function 

4.1.  Conditional  Petrophysical  Distribution  From  the 
Collocated  Data 

[is]  Both  SA-based  and  Bayesian  sequential  simulation 
approaches  require  that  GPR  velocity  be  linked  to  porosity 
through  a  conditional  petrophysical  distribution.  Although  a 
simple  approach  for  hydrogeophysical  applications  is  to 
relate  parameters  using  laboratory-derived  petrophysical 
relationships,  such  relationships  may  only  be  valid  at  the 
small  scale  and  can  encounter  problems  when  they  are  used 
to  “convert”  geophysical  images  to  subsurface  hydrologic 
properties.  Here  we  prefer  to  build  a  site-specific  conditional 
petrophysical  distribution  of  porosity  given  geophysical  data 
using  collocated  data  sets  at  boreholes.  Clearly,  in  this  case, 
the  quality  of  the  estimated  conditional  petrophysical  distri¬ 
bution  is  dependent  upon  the  number  and  quality  of  collo¬ 
cated  data  [e.g.,  Ezzedine  et  al.,  1999;  Chen  et  al.,  2001; 
Bachrach,  2006;  Paasche  et  al.,  2006].  Although  this 
approach  does  not  address  the  fact  that  the  velocity-porosity 
relation  may  in  reality  vary  in  space  due  to  differences  in  to¬ 
mographic  resolution  throughout  the  image  plane  [e.g., 
Day-Lewis  and  Lane,  2004;  Moysey  et  al.,  2005],  it  avoids 
significant  complications  and  ambiguities  associated  with 
accurately  estimating  the  spatially  variable  relationship,  and 
has  been  shown  to  be  a  reasonable  approach  to  predicting 
porosity  in  saturated  alluvial  aquifers  [. Dafflon  et  al.,  2009a]. 

[19]  In  this  study,  we  estimate  the  site  specific  condi¬ 
tional  petrophysical  distribution  between  collocated  GPR 


velocity  and  porosity  data  with  a  parametric  approach  [e.g., 
Corbeanu  et  al.,  2002;  Dafflon  et  al.,  2009a],  This  choice 
is  supported  by  the  relatively  strong  negative  correlation 
between  velocity  and  porosity,  the  somewhat  limited  size 
of  the  data  set  for  nonparametric  approaches,  and  the  exist¬ 
ing  homoscedasticity  in  the  conditional  petrophysical  dis¬ 
tribution  (Figure  3).  Based  on  preliminary  tests,  we  assume 
a  linear  relationship  with  constant  variance  between  GPR 
velocity  and  the  natural  logarithm  of  porosity;  that  is  we 
assume  that  the  conditional  expectation  of  the  natural  loga¬ 
rithm  of  porosity  E[\n{<p)\v)  given  the  velocity  v  is 
obtained  through  linear  regression  by  minimizing  the  mean 
square  error  of  the  prediction.  The  variance  cr2  (In  (<^>)  |z^)  is 
then  obtained  from  the  distribution  of  the  residuals  between 
predicted  and  measured  natural  logarithms  of  porosity. 
Based  on  Figure  3,  we  assume  that  the  conditional  petro¬ 
physical  distribution  is  approximately  Gaussian,  and  there¬ 
fore  we  only  need  to  know  the  mean  and  variance  to  draw 
values  from  the  distribution.  In  this  regard  the  natural  loga¬ 
rithm  of  porosity  is  considered  to  have  a  more-Gaussian 
distribution  of  residuals,  and  also  has  the  advantage  of 
avoiding  negative  porosity  values  in  the  simulation  process. 
Furthermore,  although  uncertainty  in  the  GPR  velocity 
tomogram  is  not  directly  included  in  this  approach,  its 
effect  is  considered  indirectly  by  the  use  of  a  field-based 
(i.e.,  in  situ)  conditional  petrophysical  distribution — which 
implies  that  the  full  range  of  porosity  values  related  to  each 
smoothed  velocity  value  is  influencing  the  petrophysical 
relationship  and  distribution  (Figure  3). 

4.2.  Spatial  Correlation  Function 

[20]  The  SA-  and  Bayesian  sequential  simulation  proc¬ 
esses  both  require  a  geostatistical  model  in  the  x-,  y-  and 


Figure  3.  Scatterplot  of  GPR  velocity  versus  the  natural 
logarithm  of  porosity  at  wells  Al,  Bl,  B2,  B3,  B4,  B6,  and 
Cl.  The  central  lines  denote  the  conditional  expectation 
E(ln(</>)|v)  assuming  a  linear  relationship  and  minimizing 
the  squared  error  when  using:  all  the  wells  (black);  all 
wells  except  B2  (red);  and  all  wells  except  Al  (green).  The 
peripheral  lines  on  the  sides  indicate  two  times  the  condi¬ 
tional  standard  deviation  cr(ln(0)|v). 
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z-directions  to  enforce  structural  variability  and  continuity 
throughout  the  realization.  In  this  study  we  estimated  the 
spatial  correlation  function  used  for  the  geostatistical 
model  by  fitting  an  exponential  model  [e.g.,  Goovaerts, 
1997 ;  Deutsch  and  Journel,  1998]  to  the  experimental  spa¬ 
tial  correlation  function  calculated  from  the  porosity  log 
data.  The  geostatistical  model  is  characterized  by  a  vertical 
range  of  4  m  based  on  porosity  log  data  from  the  full  aqui¬ 
fer  thickness ;  this  value  is  similar  to  the  vertical  range  esti¬ 
mated  for  the  full  aquifer  thickness  at  the  BHRS  previously 
by  another  approach  [Barrash  and  Clemo,  2002]  and  does 
not  produce  significant  variations  in  the  simulations  when 
changed  within  a  realistic  range  of  uncertainty. 

[21]  The  horizontal  spatial  correlation  function  was 
inferred  from  the  porosity  logs  assuming  an  isotropic  model 
because  there  is  little  evidence  for  horizontal  anisotropy. 
The  uncertainty  is  clearly  higher  in  the  horizontal  than  ver¬ 
tical  direction  because  of  more  limited  data  in  the  horizon¬ 
tal  direction.  From  the  porosity  logs  we  estimated  a  realistic 
interval  for  the  range  of  the  lateral  spatial  correlation  func¬ 
tion  to  be  between  12  and  24  m,  which  is  consistent  with 
findings  of  other  studies  for  lateral  geostatistics  of  site-wide 
hydrologic  properties  at  the  BHRS  [Barrash  and  Clemo, 
2002;  Cardiff  et  ah,  2011].  Although  the  uncertainty  of  the 
horizontal  range  could  be  included  directly  in  the  simulation 
process,  here  we  only  present  results  with  a  range  of  24  m 
because  the  findings  of  this  study  were  not  changed  much 
when  the  range  was  varied  between  12  and  24  m. 

5.  Results 

5.1.  Impact  of  Including  GPR  Data  in  Porosity 
Simulations 

5.1.1.  Predicting  Porosity  Along  GPR  Profiles 

[22]  To  assess  the  benefit  of  using  crosshole  GPR  data  to 
supplement  well  log  data  to  predict  porosity  at  locations 
where  only  crosshole  GPR  data  are  available,  a  porosity 
log  is  intentionally  withheld  from  the  simulation  process 
and  then  later  compared  with  the  prediction  at  the  withheld 
location.  This  is  done  for  two  cases:  withholding  well  B2 
and  then  Al.  The  GPR  data  consist  of  seven  profiles 
inverted  together  (Figure  1),  among  which  some  include 
well  Al  or  B2.  For  each  case,  collocated  data  for  GPR  ve¬ 
locity  and  the  natural  logarithm  of  porosity  are  used  to  infer 
the  conditional  petrophysical  distribution,  which  does  not 
change  much  when  one  porosity  log  (either  well  B2  or  Al) 
is  withheld  from  the  process  (Figure  3).  Then  100  realiza¬ 
tions  of  porosity  are  generated  for  each  case  (1)  using  sup¬ 
porting  GPR  data  with  SA-based  (Figures  4a  and  4d)  and 
Bayesian  sequential  simulation  (Figures  4b  and  4e) 
approaches,  and  (2)  using  porosity  log  data  only  with  a  se¬ 
quential  simulation  approach  (Figures  4c  and  4f ).  For  each 
case  we  also  show  the  absolute  residual  and  the  distribution 
of  the  correlation  coefficient  and  of  the  mean  absolute 
residual  between  each  simulated  and  measured  “true”  po¬ 
rosity  log. 

[23]  Use  of  GPR  data  with  either  the  SA-based  (Figure 
4a)  or  the  Bayesian  sequential  (Figure  4b)  approach  signifi¬ 
cantly  improves  the  prediction  of  porosity  at  well  B2  com¬ 
pared  to  using  only  the  porosity  log  data  (Figure  4c)  as 
seen  by  (1)  lower  residual  and  relative  errors  between 
simulated  and  true  porosity,  (2)  lower  uncertainty  in  the 


prediction,  and  (3)  higher  correlation  coefficient  between 
simulated  and  “true”  porosity.  In  particular,  by  including 
GPR  data,  the  mean  correlation  coefficient  between  pre¬ 
dicted  and  “true”  porosity  logs  increased  from  0.285  to 
0.710  with  SA-based  simulation  and  from  0.285  to  0.515 
with  Bayesian  sequential  simulation.  The  only  places  where 
GPR  data  did  not  clearly  improve  predictions  of  porosity 
are  at  the  tops  and  bottoms  of  the  tomograms,  which  are 
well  known  to  have  the  largest  uncertainty  and,  commonly, 
presence  of  artifacts  due  to  more  limited  ray  coverage. 

[24]  Prediction  of  the  withheld  porosity  log  at  Al 
(Figures  4d^lf)  also  is  significantly  improved  when  GPR 
data  are  included  in  the  simulation,  although  less  than 
at  B2  in  terms  of  the  residuals  between  simulated  and 
measured  data.  This  is  particularly  so  for  depth  intervals 
( 1)  between  838  m  and  840.2  m  where  Unit  2b  is  known  to 
have  anomalous  petrophysical  behavior  [Mwenifumbo 
et  al.,  2009],  and  (2)  between  843  m  and  844.5  m  where 
relatively  high-porosity  lenses  occur  in  Unit  4.  However, 
even  though  GPR  data  do  not  accurately  image  all  varia¬ 
tions  at  Al  compared  to  B2  (Figure  2),  still  inclusion  of 
crosshole  GPR  data  in  the  simulations  clearly  improves  the 
prediction  of  porosity  at  Al  as  demonstrated  by  nearly  dou¬ 
bling  the  correlation  coefficients  between  simulated  and 
true  porosity  (Figures  4d-4f ). 

5.1.2.  Predicting  Porosity  Outside  Crosshole  GPR 
Profiles 

[25]  Here  we  evaluate  the  benefit  of  including  GPR  data 
for  simulation  of  the  3-D  distribution  of  porosity  at  loca¬ 
tions  where  no  data  (of  any  kind)  are  available.  Again  a  po¬ 
rosity  log  is  withheld  from  the  simulation  and  then  later 
compared  to  the  predicted  porosity  at  the  log’s  location,  but 
also  the  GPR  profiles  through  this  location  are  not  included 
in  the  simulation.  Thus,  instead  of  using  seven  GPR  profiles 
as  in  section  5.1.1.,  only  two  GPR  profiles  are  included 
(B1-B3  and  Al-Cl)  when  B2  is  withheld,  and  five  are 
included  when  Al  is  withheld  (B1-B2,  C1-B2,  B2-B4,  B1-B3 
and  B6-B2).  Figure  5  shows  the  100  simulated  porosity  logs 
and  the  measured  “true”  logs  at  wells  B2  (Figures  5a-5c) 
and  Al  (Figures  5d-5f)  obtained  by  integrating  GPR  data 
with  SA-based  (Figures  5a  and  5d)  and  Bayesian  sequential 
(Figures  5b  and  5e)  simulations,  and  by  using  only  the  po¬ 
rosity  log  with  sequential  simulation  (Figures  5c  and  5f ). 

[26]  The  use  of  GPR  data  in  either  SA-based  (Figures  5a 
and  5d)  or  Bayesian  sequential  (Figures  5b  and  5e)  simulation 
still  significantly  improves  the  prediction  of  porosity  at  wells 
B2  and  Al  compared  to  using  only  the  porosity  logs  (Figures 
5c  and  5f),  even  if  no  GPR  profile  intersects  the  locations  of 
the  porosity  logs.  The  mean  correlation  coefficient  between 
predicted  and  “true”  porosity  logs  at  B2  increased  signifi¬ 
cantly  from  0.285  to  0.616,  and  from  0.285  to  0.439  by 
including  GPR  data  in  SA-based  and  Bayesian  sequential  sim¬ 
ulations,  respectively.  This  again  indicates  the  strong  contri¬ 
bution  of  GPR  data  in  constraining  the  porosity  spatial 
distribution,  even  at  locations  where  no  data  are  present. 

5.2.  Differences  Between  SA-based  and  Bayesian 
Sequential  Simulation  Approaches 

[27]  In  this  study,  the  SA-based  approach  results  in  better 
predictions  of  porosity  than  the  Bayesian  sequential  simula¬ 
tion  approach  (Figures  4  and  5).  However  universal 
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Figure  4.  Evaluating  predictions  of  measured  porosity  values  (blue  lines)  that  are  withheld  from  the  esti¬ 
mation  process  and  located  along  crosshole  GPR  profiles:  (a-c)  at  well  B2  and  (d-f)  at  well  Al.  This  is 
done  by  generating  100  3-D  realizations  (gray  lines)  of  porosity  using  crosshole  GPR  and  porosity  log  data 
in  (a  and  d)  SA-based  and  (b  and  e)  Bayesian  sequential  simulation  approaches,  and  (c  and  f)  using  only 
the  porosity  log  data  in  a  sequential  simulation  approach.  The  black  lines  show  the  means  of  all  the  simula¬ 
tions  and  the  red  lines  show  the  95%  intervals  of  occurrence.  Frequency  distributions  of  the  correlation  coef¬ 
ficients  and  of  the  mean  absolute  residuals  between  each  simulation  and  the  measured  data  are  also  shown. 


generalizations  should  not  be  drawn  because  relative 
performance  will  depend  on  details  of  a  given  system 
including  data  sets  and  relations.  The  main  difference 
between  the  approaches  is  that  in  the  SA-based  approach 


the  simulated  porosity  values  at  the  GPR  profiles  are  only 
indirectly  influenced  by  porosity  logs  through  the  spatial 
correlation  function  used  in  the  objective  function.  With 
the  Bayesian  sequential  approach  however,  simulated 
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Figure  5.  Evaluating  predictions  of  measured  porosity  values  (blue  lines)  that  are  withheld  from  the 
estimation  process  and  not  located  along  crosshole  GPR  profiles  (produced  by  removing  profiles  inter¬ 
secting  at  withheld  locations):  (a-c)  at  well  B2  and  (d— f )  at  well  Al.  This  is  done  by  generating  100  3-D 
realizations  (gray  lines)  of  porosity  using  crosshole  GPR  and  porosity  log  data  in  (a  and  d)  SA-based  and 
(b  and  e)  Bayesian  sequential  simulation  approaches,  and  (c  and  f )  using  only  the  porosity  log  data  in  a 
sequential  simulation  approach.  The  black  lines  show  the  means  of  all  the  simulations  and  the  red  lines 
show  the  95%  intervals  of  occurrence.  Frequency  distributions  of  the  correlation  coefficients  and  of  the 
mean  absolute  residuals  between  each  simulation  and  the  measured  data  are  also  shown. 
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porosity  values  are  directly  influenced  by  the  porosity  logs 
used  to  generate  the  prior  distribution  through  kriging. 
Doing  this  tends  to  create  some  averaging  between  the 
GPR  and  well  log  information,  and  thus  the  Bayesian  simu¬ 
lated  porosity  field  (e.g.,  Figures  4b  and  4e)  can  be  consid¬ 
ered  as  an  intermediate  state  between  the  SA-based  (e.g., 
Figures  4a  and  4d)  and  the  sequential  (e.g.,  Figures  4c  and 
4f )  simulation  results.  This  means  the  Bayesian  sequential 
simulation  approach  may  give  better  results  where  the 
“kriged”  prior  distribution  from  porosity  logs  (and  already- 
simulated  data)  is  reliable  and/or  where  GPR  information  is 
less  informative  (for  whatever  reasons). 

[28]  A  second  difference  between  the  approaches  is  that 
the  spatial  correlation  function  is  applied  indirectly  in  the 
SA-based  approach  through  a  global  objective  function,  but 
more  directly  influences  each  new  simulated  porosity  value 
in  the  Bayesian  approach.  However  we  did  not  observe  sig¬ 
nificant  structural  differences  in  the  results  due  to  this  differ¬ 
ence  in  how  the  spatial  correlation  function  is  applied. 
Indeed,  we  also  note  that  the  experimental  spatial  correlation 
functions  all  show  good  agreement  with  the  target  spatial 
correlation  function,  with  (Figures  4a  and  5a)  some  fluctua¬ 
tion  around  it,  and  (Figures  4b  and  5b)  slightly  smaller  mean 
of  the  spatial  correlation  ranges. 

5.3.  Accounting  For  Variable  Spatial  Statistics  and/or 
Conditional  Petrophysical  Distributions 

[29]  One  issue  about  including  GPR  data  in  the  simula¬ 
tion  of  the  porosity  spatial  distribution  is  that  the  range  of 
the  simulated  porosity  values  does  not  always  contain  the 
true  porosity  log,  as  for  example  at  A1  (Figures  4d-4f  and 
Figures  5d-5f).  In  this  regard  and  in  recognition  of  differ¬ 
ences  associated  with  Unit  2b,  which  occurs  in  most  wells 
within  about  838  and  840.2  m,  we  first  simulated  the  poros¬ 
ity  spatial  distribution  independently  for  three  layers,  with 


boundaries  at  these  depths  The  resulting  simulations  show 
that  the  porosity  predictions  were  slightly  improved  in  this 
depth  interval  at  well  A1  but  were  decreased  in  the  same 
depth  interval  at  well  B2.  In  fact  some  lateral  thickness 
and  porosity  variations  occur  between  Units  2b  and  3 
[Mwenifumbo  et  al. ,  2009]  which,  in  places,  are  not  consist¬ 
ent  with  thickness  of  the  middle  model  layer  or  are  not 
consistently  imaged  by  the  GPR  tomograms  (Figure  2). 

[30]  A  follow-up  reconnaissance  effort  was  made  to 
increase  the  variability  between  simulations  by  considering 
the  uncertainty  related  to  wells  used  to  infer  the  conditional 
petrophysical  distribution.  This  effort  consists  of  ( 1 )  calcu¬ 
lating,  independently  at  each  well,  the  likelihood  distribu¬ 
tion  for  a  range  of  parameter  values  describing  a  linear 
relationship  between  collocated  GPR  velocity  and  natural 
logarithm  of  porosity  by  grid  search  [e.g.,  Dafflon  et  al., 
2010],  (2)  summing  the  volume-normalized  likelihood  dis¬ 
tributions  obtained  at  the  different  wells,  and  (3)  drawing  a 
linear  relationship  (from  the  resulting  likelihood  distribu¬ 
tion)  to  define  the  conditional  petrophysical  distribution 
prior  to  each  new  3-D  simulation  of  the  porosity  field.  For 
this  analysis,  we  assumed  ( 1 )  the  uncertainty  in  the  meas¬ 
ured  data  to  have  a  Gaussian  distribution  with  a  standard 
deviation  of  about  0.02  and  (2)  the  constant  variance  of  the 
conditional  distribution  to  be  equal  to  the  maximum  var¬ 
iance  of  the  residual  observed  throughout  all  the  best  fitting 
linear  relationships. 

[31]  Figures  6a  and  6b  show  the  likelihood  functions 
obtained  for  a  large  range  of  evaluated  slopes  and  means  of 
natural  logarithm  of  porosity  that  define  linear  relationships, 
while  withholding  porosity  data  at  well  B2  and  then  Al, 
respectively.  Note  that  evaluating  the  mean  of  natural  loga¬ 
rithm  of  porosity  or  the  intercept  is  equivalent  for  a  constant 
mean  value  of  GPR  velocity  (here  equal  to  87.07  m/qs). 
From  Figure  6  it  is  clear  that  the  “best”  fitting  relation 


a) 


o 

CL 


cn 

o 


-1.3 

-1.4 

-1.5 

-1.6 

-1.7 

-1.8 

-0.1  -0.05  0 


slope  of  log  of  porosity  vs.  GPR  velocity 


0.2  0.4  0.6  0.8  1 

Likelihood  distribution  [] 


b)  -1.3 


c 

rn 

Ol 

E  -1.7 


-1.8 


-0.1  -0.05  0 

slope  of  log  of  porosity  vs.  GPR  velocity 


0.2  0.4  0.6  0.8  1 

Likelihood  distribution  [] 


Figure  6.  Final  likelihood  distributions  of  the  slope  and  mean  of  the  natural  logarithm  of  porosity  defin¬ 
ing  the  linear  relationship  between  GPR  velocity  and  the  natural  logarithm  of  porosity,  and  obtained  from 
independent  grid  searches  through  the  mean  square  error  for  each  of  the  wells  Bl,  B3,  B4,  B6  and  Cl,  and 
including  (a)  Al  and  (b)  B2.  Black  stars  indicate  randomly  selected  values  in  the  distribution  that  allow  def¬ 
inition  of  100  conditional  petrophysical  relationships  between  GPR  velocity  and  porosity  when  withholding 
(Figure  6a)  B2  and  (Figure  6b)  Al,  and  that  were  used  to  perform  the  simulations  shown  in  Figure  7. 
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between  GPR  velocity  and  porosity  data  changes  from  well 
to  well  (sometimes  significantly)  and  thus  the  shape  of  the 
final  likelihood  function  is  related  to  the  particular  wells 
used  in  the  analysis.  Figure  7  shows  porosity  at  B2  and 
then  A1  from  SA-based  simulation  using  GPR  data,  and 
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now  including  the  random  selection  of  a  conditional  petro¬ 
physical  distribution  from  the  likelihood  function  shown  in 
Figure  6.  Compared  to  the  previous  approach  (Figures  4a 
and  4d),  doing  this  allows  us  to  increase  the  variability 
between  the  simulations  and  thus  more  completely  include 
the  “true”  porosity  log  data.  Flowever,  while  doing  so,  this 
approach  does  not  decrease  the  mean  absolute  residual 
between  predicted  and  measured  porosity. 

5.4.  Information  Obtained  From  the  Porosity 
Simulations 

[32]  To  evaluate  prediction  of  the  overall  3-D  porosity 
distribution  at  the  BF1RS,  we  ran  100  SA-based  simulations 
using  all  31  crosshole  GPR  data  sets  together  and  all  the 
collocated  porosity  log  data.  Figure  8  shows  a  3-D  view  of 
one  randomly  selected  simulation  at  profiles  B6-B1,  B1-B2, 
B2-B3  and  B3-C3  with  horizontal  slices  at  several  depths. 
Figure  9  shows  three  randomly  chosen  realizations  at  profile 
A1-B2,  and  the  mean  and  standard  deviation  of  the  100  sim¬ 
ulations.  The  simulations  appear  realistic  and  provide  a  reli¬ 
able  3-D  estimation  of  the  porosity  distribution  (Figure  8); 
comparison  of  Figure  9  to  Figure  2  shows  that  the  relatively 
large-scale  structures  in  the  GPR  data  have  been  reproduced 
in  the  conditional  simulation,  as  has  the  small-scale  vari¬ 
ability  through  the  spatial  correlation  functions  and  porosity 
log  data.  Furthermore  the  variability  in  the  simulated  poros¬ 
ity  values  is  consistently  small  at  the  wells  and  increases 
with  distance  from  the  wells,  and  more  so  where  not  along  a 
GPR  profile.  Also,  continuity  is  well  constrained  between 
locations  that  are  conditioned  differently  (i.e.,  either  with 
well  log  data,  GPR  profiles,  or  without  conditioning  data). 

[33]  Finally,  the  3-D  simulations  of  porosity  using  GPR 
velocity  tomograms  and  porosity  log  data  (Figures  8  and  9) 
provide  very  clear  imaging  of  the:  (1)  sand  channel  (Unit 
5)  thickening  toward  the  river  at  the  top  of  the  simulation, 
(2)  relatively  high  variability  in  values  and  shorter  lateral 
continuity  in  Units  2a  and  4,  (3)  relatively  high  lateral  con¬ 
tinuity  and  lower  variability  in  Unit  3,  and  (4)  boundaries 
between  the  main  units  at  the  BPLRS.  These  imaged  details 
are  in  agreement  with  previous  geostatistical  analysis  of  the 
porosity  log  data  at  this  site  [Barrash  and  Clemo,  2002]. 

6.  Discussion  and  Conclusions 

[34]  The  main  objectives  of  this  study  were  to  develop 
modeling  approaches  to  generate  high-resolution  3-D  simu¬ 
lations  of  hydrologic  property  (porosity)  distributions  from 


Figure  7.  Evaluating  predictions  of  measured  porosity 
values  (blue  lines)  that  are  withheld  from  the  estimation 
process  and  located  along  crosshole  GPR  profiles:  (a)  at 
well  B2  and  (b)  at  well  Al.  This  is  done  by  generating  100 
3-D  realizations  (gray  lines)  of  porosity  using  crosshole 
GPR  and  porosity  log  data  in  an  SA-based  approach  and 
using  randomly  chosen  conditional  petrophysical  relations 
by  considering  the  uncertainty  in  the  relation  between  GPR 
velocity  and  porosity  (Figure  6).  The  black  lines  show 
the  means  of  all  the  simulations  and  the  red  lines  show  the 
95%  intervals  of  occurrence.  Frequency  distributions  of  the 
correlation  coefficients  and  of  the  mean  absolute  residuals 
between  each  simulation  and  the  measured  data  are  also 
shown. 
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Figure  8.  One  3-D  realization  randomly  chosen  from 
100,  obtained  from  porosity  logs  and  GPR  velocity  tomo¬ 
grams  using  an  SA-based  approach,  and  illustrated  here  at 
six  depths  in  the  aquifer  and  along  profiles  B6-B1,  B1-B2, 
B2-B3  and  B3-C3. 

geophysical  (GPR  velocity  tomography)  and  hydrologic 
(porosity  log)  data,  and  to  evaluate  how  much  benefit 
results  from  using  GPR  tomograms  with  porosity  log  data 
in  such  simulations.  Hence  we  (1)  adapted  SA-based  and 
Bayesian  sequential  simulation  approaches  to  3-D  and  for 
variable  spatial  data  coverage,  and  (2)  investigated  the 


ability  of  these  approaches  to  predict  porosity  log  data  that 
were  withheld  from  the  estimation  process  (whether  located 
along  GPR  crosshole  profiles  or  not). 

[35]  Results  show  that  use  of  GPR  velocity  tomograms  in 
the  simulation  process  improved  the  prediction  of  porosity 
compared  to  using  porosity  log  data  only.  In  particular,  by 
doing  this  (1)  the  fit  between  predicted  and  “true”  porosity 
is  improved,  (2)  the  range  of  uncertainty  estimated  by  the 
prediction  is  reduced,  and  importantly  (3)  the  correlation 
between  the  estimated  and  “true”  porosity  is  increased.  Fur¬ 
thermore  these  improvements  are  not  only  observed  when 
the  predicted  porosity  log  is  located  along  a  crosshole  GPR 
profile  but  also  when  no  supplemental  data  are  present  at  a 
given  prediction  location.  Finally,  variations  within  reasona¬ 
ble  ranges  for  parameters  involved  in  the  simulations  (such 
as  spatial  correlation  function  and  conditional  petrophysical 
distribution)  influence  local  details  of  results  to  some 
degree  but  do  not  change  the  above  findings.  For  example 
using  a  target  lateral  spatial  correlation  range  of  12  m 
instead  of  24  m  (as  was  used  here)  leads  to  (1)  small 
decreases  in  all  the  correlation  coefficients  between  the 
withheld  and  the  simulated  porosity  logs,  (2)  slightly  more 
variability  between  the  simulations,  and  (3)  an  increase  in 
the  benefit  of  using  GPR  data. 

[36]  Ultimately,  the  degree  to  which  GPR  data  can 
improve  the  prediction  of  porosity  is  dependent  on  how 
closely  the  geophysically  imaged  variations  are  related  to 
the  porosity  variations.  This  can  be  seen  for  example  with 
the  greater  improvement  in  porosity  prediction  at  B2  than 
Al,  where  the  correlation  coefficient  between  GPR  velocity 
and  porosity  is  greater  at  B2  than  Al.  This  is  not  due  to 
more  tomograms  intersecting  at  well  B2  than  Al  but  to  the 
less  successful  GPR  imaging  of  porosity  structures  at  Al. 
Thus  we  note  that  the  improvement  in  porosity  estimation 
brought  by  including  crosshole  GPR  data  can  be  reduced 
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Figure  9.  Porosity  spatial  distributions  along  profile  A1-B2,  extracted  from  3-D  realizations  obtained 
from  porosity  logs  and  GPR  velocity  tomograms  using  an  SA-based  approach:  (a-c)  three  randomly 
chosen  realizations,  (d)  mean  and  (e)  standard  deviation  of  100  porosity  simulations. 
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due  to  (1)  limited  sensitivity  to  variations  in  and  between 
some  hydrologic  units  because  of  petrophysical  behavior, 
(2)  limited  resolution  depending  on  the  acquisition  geometry 
and  physics  of  the  GPR  method,  and/or  (3)  artifacts  in  the 
GPR  tomogram.  This  means  that  the  geophysical  data  may 
be  devoid  of  value,  or  even  possibly  counter-productive,  in 
some  cases.  Possible  ways  to  improve  the  results  in  such 
cases  include  ( 1 )  consideration  of  the  spatially  variable  na¬ 
ture  of  the  relation  between  porosity  and  geophysical  data, 
(2)  use  of  supplementary  sources  of  information  such  as 
petrophysics,  amplitude  tomography,  or  electrical  tomogra¬ 
phy  to  perhaps  improve  the  characterization  of  units  where 
GPR  velocity  tomography  failed,  and  (3)  future  develop¬ 
ments  in  the  inversion  of  multiple  GPR  profiles  such  as 
with  full  waveform  inversion. 

[37]  The  SA-based  and  Bayesian  sequential  simulations 
provided  quite  similar  results  and  so  either  may  perform 
well  depending  on  the  objectives  and  requirements  of  a 
given  study.  The  main  difference  between  the  methods  as 
used  here  is  that  Bayesian  sequential  simulation  gives  more 
weight  to  the  porosity  log  data,  which  can  either  improve 
or  weaken  the  result  compared  to  the  SA-based  approach. 
Also  Bayesian  sequential  simulation  is  faster  while  the  SA- 
based  approach  allows  more  flexibility  in  accounting  for 
different  constraints  through  the  use  of  an  objective  func¬ 
tion.  Importantly  both  show  significant  flexibility  in  han¬ 
dling  the  crucial  step  of  linking  the  geophysical  data  to  the 
hydrologic  property.  And,  we  did  not  observe  significant 
changes  in  results  due  to  the  difference  in  how  the  spatial 
correlation  function  is  applied  in  the  two  approaches. 
Indeed,  both  approaches  generate  realizations  where  the  ex¬ 
perimental  spatial  functions  reproduce  well  the  parametric 
ones  used  in  the  simulation  process. 

[38]  Finally,  this  study  provides  reliable  high-resolution 
simulation  of  the  3-D  distribution  of  porosity  in  a  heteroge¬ 
neous  aquifer  at  a  real  field  site,  the  BFIRS.  As  such,  these 
results  can  be  used  further  in  examining  hydrologic  charac¬ 
terization  and  hydrogeophysical  relations  to  support  hydro- 
logic  applications.  A  next  step  will  be  to  implement  such 
porosity  reconstruction  approaches  in  more  global  strat¬ 
egies  to  predict  flow  and  solute  transport,  and  to  evaluate 
the  benefits  of  doing  so.  This  should  also  allow  improve¬ 
ment  of  our  understanding  of  the  relation  between  various 
geophysical  and/or  hydrologic  properties  in  heterogeneous 
aquifers  at  the  field  scale. 
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